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1.  INTRODUCTION 


Solving  the  equations  governing  inviscid  fluid  mechanics  is  not 
an  easy  task  - essentially  because  the  system  is  non-linear.  In  addition, 
in  steady  flow  past  blunt  bodies  various  regions  differ  from  each  other 
mathematically  - the  subsonic  flow  at  the  front  of  the  body  is  governed 
by  elliptic  partial  differential  equations  while  the  same  set  of 
equations  become  hyperbolic  farther  downstream. 


Because  of  the  difficulty  in  obtaining  analytic  solutions  there  were 
developed,  in  the  1950's,  a number  of  numerical  methods.  Two  of  the 
better  known  ones  were  the  method  of  Integral  Relations  due  to 
Dorodnitsyn  and  Belotserkovskii  [1]  and  the  Inverse  Body  Method  of 
Garabedian  [2].  Both  are  usually  used  in  the  subsonic  region  only 
with  the  method  of  characteristics  being  employed  for  the  supersonic 
flow.  These  two  algorithms  are  efficient  from  the  point  of  view  of 
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speed  of  calculation  and  core  memory  requirements.  They  do,  however, 
encounter  difficulties  in  two  areas:  near  the  sonic  point  and  at  the 

expansion  corner  (see  Fig.  1).  Since  these  methods  solve  the  steady 
flow  equations  they  cannot  be  applied  to  truly  time  dependent  problems 
such  as  the  diffraction  of  a shock  wave  by  the  bow  wave  of  a body  in 
supersonic  flight.  In  the  1960's  Lax  and  Wendroff  [3],  Richtmyer  [4] 
and  others  developed  finite-differences  algorithms  of  second  order 
accuracy  for  solving  the  time  dependent  equations.  In  principle  these 
methods  possess  several  advantages:  the  ability  to  treat  time  dependent 
problems,  the  ability  to  include  shock  waves  without  special  treatment 
and  the  fact  that  the  whole  flow  field  is  governed  by  hyperbolic 
partial  differential  equations. 

It  is  also  found  that  the  sonic  line  region  and  the  expansion 
corner  do  not  present  any  difficulty  to  these  types  of  computations. 

The  major  disadvantage  of  these  algorithms  is  having  an  additional 
dimension  (time)  - thereby  increasing  the  computation  time.  The 
length  of  computation  depends  on  the  time  step,  At,  which  the 
algorithm  allows  without  causing  numerical  instabilities.  Thus,  where 
possible,  it  is  desirable  to  devise  algorithms  with  larger  allowable 
time  step.  Zwas  [5]  has  modified  the  Richtmyer  two  step  method  so 
that  At  is  increased  by  40%  in  two-dimensional  calculations  and  by 
70%  in  three  dimensions.  Flows  containing  shock  waves  are  subject  to 
little  understood  non-linear  numerical  instabilities.  Harten  and 
Zwas  [6]  show  how  to  deal  with  this  problem  by  employing  the  Shuman 
filter.  Goldberg,  Gottlieb,  Turkel  and  Abarbanel  [7],  [8],  [9],  [10], 
developed  a number  of  algorithms  for  achieving  high  order  accuracy 
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(4th  order  and  more)  and  also  considered  the  theoretical  problems 
connected  with  applying  boundary  conditions  at  moving  boundaries  (e.g. 
the  Rankin-Hugoniot  conditions  at  a shock  wave). 

The  specific  problems  solved  in  this  report  are  as  follows: 

1.  Supersonic  flow  of  ideal  gas  past  two  dimensional  bodies  at  zero 
angle  of  attack.  The  bodies  are  blunted  wedges  connected  to 
straight  afterbodies  (see  Fig.  1).  The  calculations  were  carried 
out  for  a range  of  Mach  numbers,  2 s M £ 4,  and  for  various 
values  of  the  wedge  angle. 

2.  Supersonic  flow  past  two  dimensional  bodies  at  various  angles 
of  attack. 

3.  Supersonic  flow  past  blunted  bodies  of  revolution  such  as  blunted 
cone  followed  by  a circular  cylinder  afterbody  (see  Fig.  1). 

4.  A 3-D  calculation  of  the  flow  past  an  axisymmetric  body.  While 
this  problem  is  not  truly  axisymmetric  the  calculation  was  carried 

in  3-D  and  the  results  compared  well  with  the  axisymmetric  computations 
for  the  same  body.  [These  results  encourage  us  to  attempt  truly  3-D 
problems . ] 

In  Section  2 are  presented  the  partial  differential  equations  for 
the  various  cases;  Section  3 describes  the  numerical  scheme;  the 
boundary  and  initial  conditions  treatment  is  given  in  Section  4 and 
Section  5 discusses  the  numerical  results. 
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2.  THE  EQUATIONS  OF  MOTION 
(i)  Two  Dimensional  Flow 

The  Euler  equations  for  the  time  dependent  flow  of  inviscid, 
compressible  fluid  are: 


ft  * 37t0U)  * 37ttV>  * 0 

(Continuity) 

(2.1) 

^r(pu)  ♦ -^-(pu2+p)  + ^-(puv)  = 0 

(x -momentum) 

(2.2) 

3Y(pv)  + ^(pvu)  + gy(pv2+p)  = 0 

(y-momentum) 

(2.3) 

H + + -^rtvCE+p)]  = o 

(Energy  equation) 

(2.4) 

where  u,  v,  p,  p,  E are,  respectively,  the  fluid  velocity  in  the 
x-direction,  fluid  velocity  in  the  y-direction,  the  density,  pressure 
and  total  energy  (internal  plus  kinetic  energies)  per  unit  volume  of 
the  fluid  at  the  point  (x,  y)  at  time  t.  We  still  have  to 
characterize  the  fluid  through  its  equation  of  state.  We'll  consider 
ideal  gases  for  which 

E = + jp(u2+  v2)  (2.5) 

or,  solving  for  p, 

P - Cy-ID [E  - yf>(u2*  V2)]  (2.6) 
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where  y = cp/cv 


is  the  ratio  of  specific  heats  at  constant  pressure 
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and  volume  respectively.  For  hypersonic  flow,  for  example,  one  would 
have  to  use  a different  equation  of  state. 


The  above  system  of  partial  differential  equations  is  written  in 
divergenceless  form.  In  vector  notation  it  may  be  written  as 


3W  _ 3F  + 3 G 


3 1 3x  3y 

with  the  vectors  W,  F and  G given  by 


(2.7) 


(2.8) 


W 


p 

- m 

- n 

m 

(t-3)  m2  . (y-i)(e-^2) 

nm 

2p  2p 

p 

; f = 

; g = 

n 

ran 

4^-  n2  - (y-lHE-1"2) 

P 

2p  2p 

E _ 

m(m2.  n2)  - 

2-jJ.  n(m2+n2)  - 

_2»2  p 

2p 2 p 

and  where  m * pu  and  n = pv, 


(ii)  Three  Dimensional  Flow 

The  conservation  equations  are  of  the  form: 


3W  = 3F_  + 3_G  + 3H 
3 1 ” 3x  3y  3 z 


(2.9) 





where 


and  w is  the  velocity  component  in  the  z direction. 
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( i i i ) Cy 1 indrical ly  Symmetric  Flow 

We  get  this  case  from  the  3-D  equations  (2.9)  - (2.13)  through 
the  substitution 

x = x,  y = r sin  9,  z = r cos  e,  (2.14) 

and  by  going  from  the  velocity  components  u,  v,  w to  velocity 

components  in  the  x,  r,  9 directions  - i.e.  u , u and  u 

x r 6 

respectively.  Since  we  assume  cylindrical  symmetry,  we  may  as  well 
assume  uQ  = 0 and  label  ux  and  ur  by  u and  v respectively. 
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where  the  pressure  is  again  defined  as  in  (2.63  and  the  nonhomogeneous 
term  vector,  S,  is  given  by 


(2.20) 


The  vectors  W,  F and  G are  the  same  as  those  given  in  Eq.  (2.8). 

If  we  label  W'  £ rW  we  see  immediately  that  rF(W)  = F(W')  and 

rG(w)  - G(W').  Let  F ' = F(W')  = rF  and  G'  = G(W')  = rG;  i.e. 

F'  and  G'  are  the  same  vector  functions  of  W'  = rW  as  F and  G 

were  of  W.  Thus  our  task  becomes  the  solution  of  the  system 


3 W 1 _ 3 F ' , 3G'  Q 
3 1 3x  Dr 


(2.21) 


3.  THE  FINITE  DIFFERENCE  SCHEMES 

In  all  cases  described  herein  the  algorithms  used  are  based  on  a 
two  step  scheme  a la  Zwas  and  Burstein  [5],  [11],  [12].  In  the  two- 
dimensional  and  cylindrically  symmetric  cases  the  schemes  use  9 
computational  points  in  a 3x3  net.  In  the  3-D  case  we  require 
27  points  constituting  a 3x3x3  cube. 


We  now  present  the  schemes,  their  linear  stability  criteria  and 
the  way  we  use  the  Shuman  filter  to  prevent  non-linear  instabilities. 
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(i)  The  Two-Dimensional  Case 


The  vectors  W(x,  y,  t) , F(W(x,  y,  t))  and  G(W(x,  y,  t,))  are 

approximated  by  discretized  vectors.  Thus  Wn  = W(jAx,  kAy,  t )=W(x,y,t) 

J > k n 

where  ax,  Ay  and  Atn  are  the  step  sizes  in  the  finite  difference  net. 

Similarly  Fn  = F(W^  ) and  Gn  . = G(Wn  ).  We  shall  take 

Ay  = Ax  = constant,  but  Atn  may  vary.  The  number  of  time  steps 

n 

required  to  reach  the  time  tn  is  n ; i.e.  tn  = l Atm. 

m = l 


The  basic  finite  difference  scheme  approximating  the  system  (2.7) 
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(3.6) 


n+f 

F 

F (~(W 

+ W 

j*j,k 

2 j»i 

,k  + i 

j + i 

»k  - j 

•vn 

1 n 

n 

G 

= G(i(W 

+ W 

j + i ,k  + l 

2 j, 

. k + l 

j+l; 

. k + l 

1 n+j 

n+i 

6 

G(t(W 

+ W 

j ,k+i 

2 i+i, 

,k+i 

j -i  i 

»k+i 

(3.7) 


(3.8) 


with  similar  expressions  holding  for  discrete  vectors  with  different 
subscripts . 


The  criterion  for  the  numerical  (linear)  stability  of  the  scheme 
(3.1)  + (3.2)  constrains  the  time  step  to  be  (see  Ref.  [11]): 

AX 

At  « ■ (3.9) 

C+Xl2  + V2 

where  c = (yp/p)2  is  the  speed  of  sound.  In  practice,  one  has  to  check 
all  the  involved  qualities  at  each  grid  point  and  select  the  minimum  of 
the  right  hand  side  of  Eq . (3.9)  over  all  j and  k.  Thus,  we  use 

Ax 

Atn+1  = Vr  tn  * ' (3.10) 

max[c"ik./(uj>k)^(vyik)^ 
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(ii)  The  Cylindrical  Symmetry  Case 


For  the  sake  of  convenience  we  shall  drop  now  the  primes  of  the 
vectors  W' , F'  and  G'  appearing  in  Eq.  (2,21)  and  the  discretized 
approximations  of  these  vectors  will  be 

w”  , = W'(jAx,  kAr,  t)  = rW(jAx,  kAr,  t ) where,  as  before  Ax  * Ar 
J n n 11 

is  the  grid  size  while  tn  = £ Atm.  As  before,  we  define 

n m=l 

F,  , = F(W'(jAx,  kAr,  t ))  = F(W"  .)  and  similarly  for  G’ . 
j , k n j , k 

The  basic  scheme  representing  Eq.  (2.21)  is  given  then  by: 


n+i 

t 

j+i,k+i 


n 


, n 


W 


= w 


n 


,.n 


\ 'V;*1  f\j*1  f\jlt  Af  A 

+ ±-[F  -F  +G  -G  ]+^  S (3.11) 

j+i,k+i  2 j + 1 ,k+i  j,k+i  j+i,k+l  j + i,k  j+i,k+i 

(first  step) 


n+1  n -\,n+i  <\,n+i  ^n+i  ^n+i  ^n+i 

W = W + A [ F -F  + G -G  ] + At  S (second  step) (3.12) 

j ,k  j , k j + i,k  j - i , k j ,k+i  j,k-i  j ,k 


where 


^n  n n n 

S = — [S  + S + S + S ] 

j +i ,k+i  4 j + 1 ,k+l  j + 1 ,k  j ,k+l  j,k 


(3.13) 


-vn  + i 1 n+i  n+i  n+i  n+i 

S - T[S  +S  +S  + S ] 


j »k 


j+i ,k+i  j + i ,k- i j-i,k+i  j-i,k-i 


(3.14) 


n+i  n+i 

S - S(W  ) . 

j+i »k+i  j+i  ,k+i 


(3.15) 
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*>±,k 


*vn 

G 


G)  ,*♦{ 


etc.  are 


The  expressions  for  Fj+1>W  » rj+i>k  • uj+*,k+l  * 
the  same  as  in  the  cartesian  case,  Eqs.  (3.5)- (3.8).  The  inhomogeneous 
vector  S does  not  effect  the  (linear)  numerical  stability  and  the 
stability  condition  remains  as  in  Eq.  (3.10). 


(iii)  The  Three  Dimensional  Case 


The  various  difference  expressions  for  this  case  are  obvious 
extensions  of  the  two-dimensional  ones.  The  basic  finite  difference 
scheme  approximating  the  system  (2.9)  is  given  by: 


n+  .n  \ 'x.n  ^n 

W = W + — [F  ‘ F 

j+i,k+i,£+±  j + i ,k+i  ,£  + z 2 j+l,k+i,£+i  j,k+i,£+i 


<\,n  /\,n 

+ G - G + H - H ] 

j+j,k+l,£+j  j+i,k,£+i  j+i,k+j,£+l  j+s,k+j,£ 


(first  step)  (3.16) 


n+1 

n ^n+i  ^n+i 

W 

• V + A [F  - F + 

G 

G 

j » k , £ 

j , k , £ j + i,k,£  j - r »k  > £ 

j ,k+i  , £ 

j ,k-j  ,£ 

^n+i 

^n+i 

(3.17) 

+ H 

H ] 

(second  step) 

j ,k,£+i 

j » k , £ - i 

where 

X = At/Ax  =>  At/Ay  = At/AZ 

(3.18) 

An  1 

w = -(W  + w + w + w 

j+i,k+i,£+i  8 j + 1 ,k+l , £ j +1 »k , £ j,k+l,£  j,k,£ 


+ W + W + W +W  ] (3.19) 

j + l,k+l,£+l  j+l,k,£  + l j,k+l,£+l  j »k , £+1 
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'x.n  i n n n n 

F = F(r(w  +w  +W  +W  )) 

j+l,k+i,£+i  j+l,k+l,£+l  j +1 , k+1 , £ j+l,k,£+l  j + 1 , k , £ 


(3.20) 


n + i 


n+i 


n + i 


n+i 


n+i 


F = F (4-(  W + W + W + W )) 

j +i * k , £ j + i ,k+i , £ + i j + i , k+i , £ - i j+i,k-i,£+i  j+i,k-i,£-i 

(3.21) 


n 


n 


n 


1 n 

G = G (— (W  + W + W + W )) 

j+i,k+l,£+i  4 j+l,k+l,£+l  j +1  ,k+l , £ j,k+l,£+l  j,k+l,£ 


(3.22) 


<\,n+i 

G 


j »k+i , £ 


i n+i 

= GCT(w 


+ W 


n+i 


n+i 


W 


n + i 


+ W 


)) 


j+i,k+i,£+i  j+i*k+i,£-i  j_i,k+i,£+i  j*i»k+i,£-j 

(3.23) 


<\.n 

H 


1 n 
= H(-(W 


n 


n 


+ W 


+ W 


+ w 


j + 5 ,k+i ,£  + 1 H j+l,k+l,£  + l j + 1 ,k  , £ + 1 j,k+l,£  + l j,k,£  + l 


)) 


(3.24) 


^n+i 

H 


n+i 


n+i 


n+i 


j >k, £ + i 


i n+i 

H(— (W  + W W + W ) ) , 

4 j +i ,k+i , £+i  j+i,k-i,£  + i j -i ,k+i , £ + i j-itk-i,£+i 


(3.25) 


with  similar  expressions  for  different  subscripts.  By  anology  to  (3.10) 
the  largest  time  step  allowable  under  the  (linear)  stability  criterion  is 


Ax 


At  . , * t_  . , - t $ 


"n+1  n+1 


(3.2.6) 
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Civ)  The  Treatment  of  Shock  Waves  by  an  Automatic  Numerical  "Switch" 

In  most  of  the  flow  field,  the  2-D  and  3-D  algorithm  - Eqs.(3.1)+ 

(3.2)  and  ( 3 . 16) + (3 . 17)  and  the  cylindrical  symmetry  algorithm  - Eqs. 

(3 . 11 ) + (3 . 12)  - give  results  which  are  linearly  stable  and  which  are  of 
second  order  accuracy.  In  the  vicinity  of  shock  waves  and  stagnation  points 
there  exists  (for  different  reasons)  the  danger  of  non-linear  numerical 
instability.  Harten  and  Zwas,  [6],  ameliorate  this  phenomenon  by  a 
modified  application  of  the  Shuman  filter.  Usually  (see  Vliegenthart  [13]) 
the  filtering  is  applied  to  the  whole  flow  field  and  this  reduces  the 
accuracy  of  the  algorithm  to  first  order.  If,  however,  the  filtering  is 
done  only  in  the  immediate  vicinity  of  a shock  wave,  then  the  non-linear 
instability  is  usually  prevented  while  the  accuracy  of  the  computation  in 
the  rest  of  the  flow  field  remains  of  second  order. 

In  the  two-dimensional  and  the  cylindr ically  symmetric  cases 
one  proceeds  as  follows: 


n+1 

n 

W 

L W 

(3. 

27) 

j .k 

j y k 

n+1  n+1  i 

X 

n+1 

n + 1 

X 

n+1 

n+1 

f = W + - 

j ,k  j ,k  4 

[e 

(W 

• W ) 

- 0 

(W 

- w )] 

j+s  >k 

j + l,k 

j >k 

j ' 

. k j ,k 

j -l.k 

1 

y 

n+1 

n+1 

y 

n+1 

n+1 

+ — 

4 

[e 

(W 

• W ) 

- 0 

(W 

- W )]  . 

(3. 

28) 

4 

j »k+i 

j ,k+l 

j »k 

j : 

>k- 

■i  j >k 

j » k- 1 

Where  the  operator  L is  the  scheme  (3.1)+(3.2),  or  in  the  cylindrical 
symmetry  case,  the  scheme  (3 . 11) + (3 . 12) . The  "switches"  0^,  0X 
are  defined  as  follows: 
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0 = x 
j + i,k 


T I p - p I 

j4l»k  j-»k 

max 

L j > k 


p - p 

j+l,k  j,k  J 


(3.29) 


j »k+i 


p - p 

j>k+!  J>k 


-*  m 


max 

j >k 


p - P 

j,k+l  j,k 


(3.30) 


Near  shock  waves,  or  other  regions  of  very  strong  gradients,  the 

expressions  in  the  square  brackets  in  (3.29)  and  (3.30)  are  of  order 

unity  and  then  0X,  ey  : x and  the  filtering  defined  by  (3.28)  becomes 

operative.  Away  from  the  shock-wave,  the  flow  is  smooth  and  ex*^  z 0[(Ax)m] 

=n+ 1 n+ 1 

Thus,  for  mil,  in  the  smooth  regions  W = W + 0(ax2)  at  least; 

j#k  j»k 

i.e.  second  order  accuracy  is  preserved.  In  practice,  one  uses  the 
scheme  (3 . 27)  + (3 . 28)  with  the  0X's  and  0^'s  substituted  from  (3.29) 
and  (3.30).  Because  of  linear  stability  requirements  we  are  constrained 
to  use  0 < x $ 1. 


n+1 


W 


n+1  i x 

w . * tK 


i y 

+ 


lr  2 
+ -[9 


case  Eq. 

, (3.28) 

takes  the 
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x y 

where  0 and  9 are  defined  respectively  according  to 

j,k+i,n 

(3.29)  and  (3.30)  suitably  modified,  and 


0 


z 

j 


r I p * p I -im 

j,k,&+_l 


max 
*-j  ,k,«, 


p - p 

j,k,t  + l j,k,£ 


(3.31) 


Also,  linear  stability  analysis  [6],  shows  that  in  contradistinction  to 
the  two  dimensional  and  cylindrical  symmetric  cases,  in  the  three 
dimensional  case  the  filtering  coefficient  x has  a more  restricted 
range.  Specifically,  in  the  3-D  case  we  are  constrained  to  use 
0 < x S 2/3. 

4.  TREATMENT  OF  BOUNDARY  AND  INITT.AL  CONDITIONS 
(i)  Boundaries  That  Are  Not  on  The  Body 

The  computation  is  usually  done  over  a rectangular  grid  of  J x k 
net  points,  where  J is  the  number  of  grid  points  in  the  x-direction  and 
K is  the  number  in  the  y or  r direction.  We  choose  K in  such  a way 
that  the  bow  shock  wave  will  not  cross  the  upper  boundary,  k - K (the 
lower  boundary,  k = 1,  is  usually  taken  along  the  axis  of  symmetry)  but 
the  right  hand  boundary,  j = J.  (See  sketch.) 


Along  BC  the  boundary  conditions  are  found  by  extrapolation  along 
45°  lines  except  that  very  near  to  B (2  points  along  BC)  where 
the  extrapolation  is  in  a direction  perpendicular  to  BC . Along  CF 
we  use  the  same  strategy  except  that  very  near  F the  extrapolation 
is  in  a direction  parallel  to  the  body  surface.  Along  AE  the 
boundary  conditions  are  determined  by  the  symmetry  of  the  flow  (zero 
angle  of  attack). 
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If  the  flow  is  not  symmetric  about  the  x-axis  (angle  of  attack  is 
not  zero)  then  for  the  lower  boundary  we  use  not  AE  but  a line  B'C' 
which  we  treat  in  the  same  manner  as  BC . On  the  boundary  AB  we  set 
fixed  the  ambient  free  stream  conditions  of  the  steady  state  flow 
which  we  are  trying  to  model. 

(‘i i ) Boundary  Conditions  on  the  Surface  of  a Two-Dimensional  Body 

One  way  of  dealing  with  boundary  conditions  on  a body  of  arbitrary 
shape  is  to  transform  the  computational  grid  in  such  a fashion  that  the 
body  surface  then  coincides  with  one  of  the  new  coordinates.  The 
difficulty  with  this  is  that  the  finite  difference  algorithm  becomes 
more  complex  and  has  to  be  changed  to  fit  each  new  problem.  We  chose,  on 
the  other  hand,  to  stay  with  the  convenience  of  a rectangular  mesh.  Of 
course,  we  then  face  the  problem  that  the  body  surface  does  not,  in 
general,  pass  through  grid  points.  (See  sketch  on  previous  page.) 

We  need  to  know  the  components  of  W at  the  point  "Q"  (inside  the 
body)  at  time  tn  in  order  to  be  able  to  compute,  for  example,  W at 
point  "2"  (outside  the  body)  at  time  tR+^  “ t + Atn*  The  points 
"a"  and  "b"  were  chosen  in  such  a way  that  the  line  Qab  is  normal  to 
the  body  at  the  point  of  intersection,  "c".  For  the  purpose  of  the 
discussion  in  this  section  only,  let  f be  any  component  of  W or  of 
the  related  vector  (p , u,  v,  E) . When  f stands  for  either  the  density, 
energy  or  the  velocity  in  the  direction  of  the  tangent  to  the  body  at  "c" 
we  find  its  value  at  Q by  using  a "parabolic  reflection".  Namely,  we 
pass  a parabola  through  the  points  a,  b and  Q so  that  the  derivative 
in  the  direction  of  the  normal  Qab  is  zero  at  the  point  "c".  This  is 
done  by  setting 


19 


- f 


_a 

2 


u; 


* ) + f , 

aJ  a ’ 


(4.1) 


where  i , d,  and  nn  are  the  distances  of  the  points  "a",  "b"  and 


’Q”  from  the  point  "c”.  The  values  of  f and  f,  are  found  from 


parabolic  interpolation  among  the  grid  points  nearest  to  them.  Thus 


fa  * y(a-l)(a-2)f1  - a (o  - 2 ) f 2 ♦ *1^11  f , 


(4.2) 


where  aAy  is  the  distance  of  point  "a"  from  point  "1",  etc. 

When  f stands  for  the  velocity  in  the  direction  normal  to  the  body, 
it  must  vanish  at  point  "c".  We  satisfy  this  condition  through  a 


parabolic  extrapolation  that  yields  for  the  value 


f f, 

fq  ' V‘i7  ‘ A 


-)  i 


& . 

a - b 


i f 
a + _ii  . 

0 J 


(4.3) 


(iii)  Boundary  Conditions  on  the  Surface  of  a Body  of  Revolution 


The  philosophy  of  the  treatment  is  the  same  as  in  Section  4-  (ii) 
except  that  where  the  radius  of  curvature  of  the  body  is  finite  we  use 
instead  of  parabolic  reflection  and  extrapolation  linear  ones.  This 
helps  with  the  stability  and  leaves  the  overall  accuracy  unchanged. 
Thus  we  replace  (4.1)+(4.2)  by  (4.4)+(4.2): 


fQ  = fa  + l 


aQ  ' £a 


C£b  - £a) 


(4.4) 
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while  (4.3)  is  replaced  by 


(4.5) 


Note  that  even  though  the  finite  difference  system  is  solved  for 
W'  = rW,  the  conditions  (4.4)  and  (4.5)  are  applied  to  the  physical 
qualities  W * W'/r. 

Near  the  axis  of  symmetry,  r = 0,  we  have  the  problem  of  W'  = 0 
there.  To  compute  W on  the  axis  we  use  the  known  values  at  r = Ar 
r = 2Ar  and  r = - Ar  and  interpolate.  Finally,  note  that  also  on 
boundaries  away  from  the  body  surface,  such  as  BC  for  example,  all 
extrapolations  are  done  on  W and  not  on  W’. 

( iv)  Initial  Conditions 

At  t = 0 the  whole  flow  field  is  assigned  the  free  stream  conditions. 

We  chose  to  nondimensionalize  in  such  a way  that  both  the  free  stream 

pressure  and  density  take  on  the  value  of  1.  Thus  the  free  stream  sound 

speed  becomes  c = fy. 

00 

When  we  did  parametric  runs  the  conditions  at  t = 0 were  set  to 
the  converged  solution  of  a similar  run  thus  saving  computation  time. 
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5.  NUMERICAL  RESULTS 

The  numerical  results  were  obtained  for  several  problems. 

(i)  Steady,  two  dimensional,  supersonic  flow  past  a circularly  blunted 
wedge  with  a semi-apex  angle  of  13°,  at  zero  angle  of  attack. 

(ii)  Steady,  cy 1 indrical ly  symmetric,  supersonic  flow  past  a spherically 
blunted  cone  with  a semi-apex  solid  angle  of  13°,  at  zero  angle  of 
attack . 

In  both  of  the  above  groupings  the  computations  were  carried  out 
for  free  stream  Mach  number  range  of  2 $ Mk  s 4 with  jumps  of 
AM^  =0.5  from  one  run  to  another.  The  graphs  show  the  distribution 
along  the  body  of  the  ratio  of  surface  pressure  to  the  stagnation  point 
pressure  and  the  distribution  of  the  surface  local  Mach  number. 

The  Mach  number  distribution  over  the  wedge  is  shown  in  Figs.  2-6. 
The  pressure  distribution  over  the  wedge  is  shown  in  Figs.  7 - 11.  The 
pressure  distribution  over  the  cone  is  shown  in  Figs.  12  - 16.  The  Mach 
number  distribution  over  the  cone  is  shown  in  Figs.  17  - 21.  At  the  top 
of  Figs.  2,  7,  12  and  17,  each  at  Mb  = 2,  is  shown  the  body  shape  over 
which  the  calculation  was  done. 

The  surface  pressure  over  the  body  was  computed  in  two  ways:  directly 
from  the  finite  difference  scheme  and  also  by  assuming  that  the  body 
represents  a stream  tube  over  which  there  is  isentropic  flow  and  hence 
the  pressure  over  it  is  related  directly  to  the  (local)  surface  Mach 
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number.  In  Fig.  12,  for  example,  we  show  the  pressure  as  computed  by 
both  approaches.  The  dashed-line  curve  gives  the  pressure  ratio  as 
obtained  directly  from  the  finite  difference  equation  and  the  undashed 
curve  corresponds  to  the  "isentropic"  calculation.  It  is  seen  that  the 
results  are  nearly  identical  except  near  the  front  of  the  cone  where 
the  calculations  are  affected  by  the  small  r value.  Because  of  the 
agreement  between  the  two  methods,  we  show  on  most  graphs  only  pressure 
distribution  curve. 


All  the  two  dimensional  calculations  were  done  on  a grid  of 
52x55  (J  = 52,  K = 55).  Running  time,  when  the  initial  conditions 
correspond  everywhere  to  the  free  stream  value  is  about  25  minutes 
(there  are  some  variations  depending  on  Mach  number,  wedge  angle,  etc.). 
But  if,  for  example,  for  the  = 2.5  run  we  use  as  initial  conditions 

the  numerical  solution  from  the  Ma  = 2.0  run,  then  the  running  time 
decrease  to  about  10  minutes.  We  thus  found  that  the  average  running 
time  per  case,  for  computing  the  cases  = 2,  2.5,  3,  3.5,  4 is 

about  12  minutes. 


For  the  flow  around  the  blunted  cone  we  used  a net  of  65x64  grid 
points  and  the  computation  time  was  roughly  the  same  as  in  the 
cartesian  case. 


In  order  to  compare  our  algorithm  with  other  numerical  techniques 
we  made  use  of  the  results  obtained  from  semi-empirical  computer  programs 
based  on  Russian  data  for  a blunt  cone  with  10°  semi-apex  angle  at 
M = 3.  This  information  is  contained  in  a 1966  AVCO  Report 
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no.  SR  10-TR-66-47  written  by  R.ll.  Kohrs.  We  ran  calculations  for  the 
same  conf igurat ion . The  comparison  is  shown  in  Fig.  22.  It  is  seen 
that  the  agreement  is  very  good. 

(iii)  In  Fig.  23  we  show  the  results  for  our  blunted  wedge  but  at  an 
angle  of  attack  of  5°.  The  distributions  are  shown  for  both  the 
upper  and  lower  surfaces. 

In  all  of  the  above  runs  we  used  a linear  Shuman  filter,  i.e.  we 
took  m = 1 in  equations  3.18  and  3.19.  The  dissipation  coefficient 
X was  taken  to  be  1/2  in  the  two  dimensional  calculations.  In  the 
axisymmetric  case  the  value  of  x was  varied  to  get  best  results  for 
the  stagnation  density  and  was  found  to  be  . 9 s x < 1.0. 

(iv)  Finally,  we  tested  our  3-D  package  by  applying  it  to  the  problem 
of  the  supersonic  flow  past  a body  composed  of  a hemisphere 
followed  by  a circular  cylinder,  at  a zero  angle  of  attack. 

This  allowed  us  to  check  how  the  results  obtained,  using  a 3-D 
algorithm  compare  with  those  given  by  a (two-dimensional)  axisymmetric 
scheme.  The  computational  net  was  42x40x40.  Thus  we  had  67,200  mesh 
points  as  compared  to  the  4,225  points  of  the  65x65  "2-D"  mesh.  In 

addition,  in  each  mesh  point  in  the  3-D  case  we  have  to  store  a 
5-vector  (p , pu,  pv,  pw,  E)  as  compared  to  the  4-vector  (p , pu,  pv,  E) 
in  the  axisymmetric  case.  Thus  the  storage  requirements  in  the  3-D  case 
exceed  by  a factor  of  20  (twenty)  those  of  the  2-D  case.  Since  this 
requirement  exceeds  the  core-memory  capacity,  we  used  discs  for  the 
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inass  - storage . Here  we  utilized  the  hyperbolic  nature  of  the  p.d.e.  system 
as  each  field  point  was  computed,  its  "cube  (3x3x3)  of  influence"  was 
moved  by  one  mesh  point  freeing  core-memory  storage  for  data  to  be 
transfered  from  the  disc.  The  data  transfer  can  be  done  while  the 
arithmetic  unit  carries  out  the  computation.  In  this  manner  the  effect 
of  the  slow  rate  of  transfer  is  mitigated.  In  fact,  a typical  run  took 
12  times  longer  than  the  corresponding  2-D  calculation  (all  with  the 
above  given  mesh  sizes).  The  "improved"  efficiency  (12  vs.  20)  is  due 
to  the  coarser  mesh  (1/40  vs.  1/65).  The  pressure  distribution  thus 
obtained  agrees  well  with  the  axisymmetric  results.  Typically,  while 
the  stagnation  pressure  was  under-predicted  by  about  3%  in  the  axi- 
symmetric calculation,  it  was  over-predicted  by  about  4%  in  the  3-D 
runs.  The  drag  coefficient 
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was  calculated  in  both  cases.  Typical  values,  at  = 3,  are 

Cp  = 0.98  (axisymmetric) 

CD  = 1.04  (3-D)  . 

We  conclude  therefore,  that  our  3-D  algorithm  is  apparently  reliable 
and  we  now  plan  to  apply  it  to  truly  three  dimensional  flow,  i.e.  in 
the  case  of  non-zero  angle  of  attack. 
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comparison  of  the  pressure  distribution  along 

A BLUNTED  CONE  AS  COMPUTED  BY  FINITE-DIFFERENCE 
METHOD  (DASHED  CURVE)  WITH  SEMI-EMPIRICAL 
FORMULA  (OPEN  CIRCLES). 
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